function [c_new,dt_required]=RK45(sourcefunc,c_old,dt,tolerance,start_time,safty_factor)

%% ODE system solver Runge-Kutta adaptive [4-5]
% Code by Kaveh Zamani
% Civil and Environmental Engineering, UC Davis
% Method of Cash, J.R., and Krap A.H. (1990) See Press book 16.2

% u_new: solution at t = n + 1
% u_old: value at t = n
% delta: desired accuracy
% error
% safty_factor = 0.95
%
%
%%


k1 = dt*feval(sourcefunc,start_time          ,c_old       );
k2 = dt*feval(sourcefunc,start_time +   dt/5 ,c_old + k1/5);  
k3 = dt*feval(sourcefunc,start_time + 3*dt/10,c_old + 3*k1/40 + 9*k2/40);
k4 = dt*feval(sourcefunc,start_time + 3*dt/5 ,c_old + 3*k1/10 - 9*k2/10 + 6*k3/5);
k5 = dt*feval(sourcefunc,start_time +   dt   ,c_old -11*k1/5  ...
                                         + 5*k2/2 - 70*k3/27 + 35*k4/27);
k6 = dt*feval(sourcefunc,start_time + 7*dt/8 ,c_old + 1631*k1/55296 + 175*k2/512 ...
                                          + 575*  k3/13824  +    ...
                                            44275*k4/110592 +    ...
                                            253*  k5/4096      );
                                        
error = abs(                               ...
            (37/378  -  2825/27648)* k1  + ...
            (250/621 - 18575/48384)* k3  + ...
            (125/594 - 13525/55296)* k4  + ...
            (        -   277/14336)* k5  + ...
            (512/1771-     1/4    )* k6  )  ;
    
delta = abs(tolerance*c_old) ;

if (error <= delta)
  dt_required = dt*safty_factor* (delta/error)^(1/5);
else
  dt_required = dt*safty_factor* (delta/error)^(1/4);  
end

c_new = c_old + 37*k1/378 + 250*k3/621 + 125*k4/594 + 512*k6/1771;

return
end